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Abstract 

We investigate the optical properties of layered structures with graphene at the inter- 
face for arbitrary linear polarization at finite temperature including full retardation by 
working in the Weyl gauge. As a special case, we obtain the full response and the related 
dielectric function of a layered structure with two interfaces. We apply our results to dis- 
cuss the longitudinal plasmon spectrum of several single and double layer devices such as 
systems with finite and zero electronic densities. We further show that a nonhomogeneous 
dielectric background can shift the relative weight of the in-phase and out-of-phase mode 
and discuss how the plasmonic mode of the upper layer can be tuned into an acoustic 
mode with specific sound velocity. 

Pacs: 78.67.Wj, 73.21.Ac, 42.25.Bs, 73.20.Mf 

1 Introduction 

Graphene, the two-dimensional allotropc of carbon, has become one of the most active 
fields in today's both experimental and theoretical condensed matter physics[l, 2, 3, 4]. 
Among many extraordinary phenomena and proposals, the optical properties of graphene 
have generated particular interest due to potential applications, [5, 7, 8, 6] but also because 
they are intimately related to the discovery of exfoliated graphene. [9] 

Recently, plasmonics based on graphene has become a new emerging subficld, trying 
to take advantage of the strong electronic confinement and long propagation lengths of its 
carriers and, most importantly, the possibility of applying an electrostatic gate voltage. [10, 
11, 12, 13, 14, 15] Due to the linear spectrum, different gate voltages can have strong 
effects on the carrier concentration and thus on the plasmonic spectrum, especially for 
small nano-islands.[17, 16] 

There is also renewed focus on layered structures due to the experimental advances 
of exfoliating a number of materials and placing them on top of each other. [18, 21, 19] 
Like this, the Coulomb drag of closely separated graphene layers was observed, [20] which 
has triggered considerable attention by various theoretical groups. [22, 23, 24, 25, 26, 27, 
28, 29, 30] In this paper, we want to combine these two fields and analyze the plasmonic 
spectrum of layered structures at finite temperature. 

To discuss the undamped plasmon dispersion, it suffices to determine the zeros of the 
dielectric functions[31] or, alternatively, the zeros of the denominator of the transmission 
or reflection amplitude. [32] But for a finite imaginary part, i.e., especially for finite tem- 
peratures, this procedure is ambiguous. We will thus need a different approach and will 



investigate the energy loss function. For single layer graphene systems, the energy loss 
function is related to the (negative) imaginary part of the inverse of the dielectric func- 
tion, but for double layer structures, this function changes its sign. [32] We will thus define 
the energy loss function as the trace of the (negative) imaginary part of the full response. 

Another common approximation to obtain the longitudinal plasmon dispersion consists 
in using the static Coulomb potential and for small wave numbers q <C kp in the regime 
of small layer separation k F d <C 1, the Coulomb potential can be further simplified, only 
depending on the average of the outer dielectric media (ei + e^)/2. For the general case, 
the full electrostatic problem has to be considered, which has been recently done in the 
context of graphene double layers. [32, 33, 34] 

Since for small energies of the order of ae F {a being the fine-structure constant) re- 
tardation effects lead to strong light-matter interaction, [35] we will here derive the full 
retarded photon propagator for the longitudinal and transverse channel. Another focus 
of this work is placed on materials with a large dielectric constant which strongly screen 
the graphene layers. SrTi03, e.g., has a relative dielectric constant of e <~ 300, which can 
reach up to ~ 5000 at liquid helium temperature due to the proximity of a ferroelectric 
instability. [36] Another example are surface states of a three-dimensional topological insu- 
lator which are separated by the width of the sample. For Bi 2 Te3, the two electronic Dirac 
systems are then electrostatically coupled through a dielectric medium with e <~ 100. [33] 
Particularly, wc will show that strong dielectrics shift the relative weight of the in-phase 
and out-of-phase mode. Used as a substrate, they strongly screen the graphene layers and 
therefore basically act like metals, leading to a linear plasmon dispersion. [37, 38] 

The paper is organized as follows. In section II, we derive the photon propagator in 
free space, separating the final result into longitudinal and transverse channels. In section 
III, we formulate the linear response theory for a layered structure including graphene in 
the interfaces and give explicit expressions for the double layer geometry. In section IV, we 
finally discuss the near-field optical properties of layered graphene structures, contrasting 
between single and double layer, high- and low temperature, large and small dielectric 
constants. We close with a summary and an outlook. An appendix outlines the analytical 
discussion how to obtain the linear plasmon dispersion for a double layer system with 
large-e substrate. 

2 Photon propagator in a homogeneous medium 

In this section and the following section, we will derive the retarded photon propagator 
in a homogeneous and nonhomogeneous medium, outlining all details. For alternative 
introductions to confined photon systems, we refer the reader to Ref. [39]. Even though 
our treatment will be entirely classical based on Maxwell's equations, we will still call 
the final result, Eq. (15), the photon propagator since this expression coincides with 
the quantum mechanical photon propagator. [40] The reader only interested in the results 
presented in Sec. IV may skip this and the following section. 

We depart from Ampere's circuital law equation including Maxwell's displacement 
current assuming a harmonic time evolution with frequency to: 

V x H(r) = -iwD(r) + j(r) (1) 

We will first recall the usual representation of the photon propagator in Cartesian coordi- 
nates and then introduce the representation in cylindrical coordinates, suitable to discuss 
layered structures. 



2 



2.1 Photon propagator in Cartesian coordinates 

Introducing the electrostatic potential </> and the vector potential A by 

E(r) = iwA(r) - V*(r) (2) 
H(r) = —V x A(r) (3) 

we obtain from Eq. (1) the following equation: 

V x V x A + kj^ e£o(i^A - V</>) = fifi j (4) 

2.1.1 Lorentz gauge 

In the Lorentz-gauge V • A = iojfifioesocf), Eq. (4) can be written as four independent 
inhomogeneous Helmholtz equations 

(-V 2 -fcg)A"(r)=i"(r) (5) 

with the four-dimensional vector potential — (<fr, A) and the four-dimensional current 
= (p/eeo, Woj)- The dispersion relation reads 

, 22 uj 2 uj 2 

where c denotes the speed of light in vacuum and c\ — c/y/JIe the (slower) speed of light 
in the dielectric medium characterized by fi and e. The Green's function defined by 

(-W 2 -k 2 )G (r,r') = 5(r-r') (7) 

is thus a scalar and reads 

p ±ifco|i — r'\ 

G (ry) = — (8) 

47rjr — r'\ 

where the plus-sign defines propagation out of the source and the minus-sign convergence 
into the source. Within the Lorentz gauge, the electric field due to a given current density 
j, as defined in Eq. (2), is thus given by 



E(r) = iw + ^ ) A(r) = i^Mo (l + ^f) J d\'G {r, r')j(r') 



(9) 



2.1.2 Weyl gauge 



In the following, we will not work in the Lorentz gauge, but will set the electrostatic 
potential equal to zero, i.e., <p — 0. This gauge condition is often refered to as the "Weyl 
gauge". The Weyl gauge implies that we only need to consider the propagation of the 
vector potential and is often used when interactions with non-relativistic particles are 
involved. Graphene's linear response to the incoming light field is thus entirely defined by 
its current density. With = 0, Eq. (4) becomes 

V x V x A - k 2 A = /Li/ioj • (10) 

In this gauge, the operator acting on A is a vector that connects different spacial directions. 
For a general solution, we thus need to determine the dyadic Green's function defined by 

(V x V x -k 2 )G = IS(r - r') , (11) 
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where I is the 3x3 unit tensor. At first glance, it seems that the Weyl gauge results in 
more complicated expressions of the Green's function. This is true in real space, but in 
Fourier space, the expressions are only slightly more involved. And by having eliminated 
the electrostatic potential </>, the boundary conditions of the subsequent scattering problem 
become more compact since they only have to be satisfied by the vector potential A. 

We can easily obtain the dyadic Green's function by noting that the Maxwell's equa- 
tions yield the following expression for the electric field: 

V x V x E - k 2 E = iwfifioj (12) 

The electric field is thus defined by the same dyadic Green's function as the vector potential 
in the Weyl gauge. Since the electric field is gauge independent, we can use the expression 
of Eq. (9) to deduce the dyadic Green's function G from the scalar Green's function G . 
In real space, we obtain 



G(r,r')=(l+^)G (r,r') (13) 



and in Fourier space 



/ k a k"\ 

G a '?(k,uj) = (6 a ,0 - -^-j G (k,u) . (14) 

With Go{k,iu) = (k 2 — fcj}) -1 , we finally obtain the retarded photon propagator T>q^ = 
—fifioG™' 13 in the Weyl gauge, 

r,<*P{u \ W fx k a k !i \ 

Decomposing it in longitudinal and transverse components, it reads 

vf(k,co) = Dl^f + D" T (<W - ^pj (16) 

with the functions DL^(k,uj) given by 

n" - 1 no _ W^o fi7l 

££ W fc - [iO/CiY 



2.2 Photon propagator in cylindrical coordinates 

In a layered structure, assumed to be perpendicular to the z axis, the components par- 
allel to the interface, q — (q x ,q y ), will be preserved as a good quantum number. It is, 
therefore, convenient to employ the following representation for the Green's function in a 
homogeneous medium: 

Vf(z,z';q,uj) = ^ J dk z e ik ^^ Vf(k, to), (18) 

with k = (q,k z ). 

In this representation, the tensor components have a different structure depending on 
whether a,(3 = i, j with i, j = x, y or a, (3 — z. For the in-plane components of the tensor 
2?q'' 3 , we have 

V«(z,z>-,q,u) = -^e-^-*'\ (V^f) (19) 
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with q' = yjq 2 — (w/ci) 2 . Decomposed into longitudinal and transverse contributions, we 
obtain 



Vi j (z,z') = d°e-^- 2 'l ^ +d° t e~ q ' lz ~ z ' 1 ^ - ^ 
with the in-plane propagators df t (q,ui) given by 



(20) 



</•/ J ':!':;■ w 

a,0 



2ee cj 2 ' * 2q' 

For the cross terms of T>q' p , we get 



T>i*{z, z') = Vl\z, z') = df-^e-^ z - z '\ sgn(z - z') . (22) 

This shows that in-plane longitudinal sources can generate not only in-plane, but also out- 
of-plane fields with a phase shift of n/2. On the other hand, in-plane transverse currents 
can only generate in-plane transverse fields. 

Finally, out-of-plane sources generate out-of-planc fields given by the tensor component 



p ( j--(-.. :')=,/}'( 2*^^ -^e-^l"-' 'I) . (23) 



3 Linear response of a layered geometry including graphene 

We now consider the experimentally important situation of layered structures, i.e., we 
assume the existence of well-defined interfaces at which the material properties are discon- 
tinuous. The Maxwell equations written in integral form then yield boundary conditions 
for the normal and tangential field components, see e.g. Ref. [41]. For the normal compo- 
nents they read 

n-(D 2 -Di)=p, n-(B 2 -Bi) = (24) 

with p the charge density on the interface. 
For the tangential components they read 

nx(E 2 -E 1 ) = 0,nx(H 2 -H 1 )=j (25) 

with j the current density on the interface. 

An arbitrarily polarized electromagnetic wave can always be expressed by superpos- 
ing two linearly polarized waves which are orthogonal to each other. We can thus define 
p-polarized and s-polarized waves, respectively, according to the plane of incidence. Al- 
ternatively, we will use the denomination of longitudinal and transverse polarization. 

The boundary conditions for the normal and tangential field components arc not in- 
dependent of each other since they are connected by Maxwell's equations. According to 
the plane of incidence, we will either use Eqs. (24) in the case of longitudinal polariza- 
tion or Eqs. (25) in the case of transverse polarization. These conditions then simplify 
considerably when working in the Weyl gauge. 

Notice that the influence and properties of graphene are only accounted for via the 
charge and current densities on the interface, p and j. The properties of graphene thus 
enter when matching the vector field at the interface of two adjacent dielectric media. 
Since we have set = 0, these properties are entirely contained in the current-current 
response of graphene, Xij- The superindex denotes the bare current response, i.e., the 
response to the total (external plus induced) vector potential. 



5 



We will now assume an isotropic system such that current-current response tensor can 
be split up into a longitudinal and a transverse contribution, i.e., 

x! j =x!f+x < l(6 ij - 9 f-). (26) 

Within the Dirac approximation this decomposition is always possible and only for transi- 
tions close to the van Hove singularity the full tensor structure needs to be considered. [45] 
The full photon propagator in the presence of a single graphene layer at the location 
z\ modifies the "vacuum" propagator in the following way: 

V^(z,z')=V^(z,z')+Vf(z,z 1 ) Xij 'Di (z 1 ,z') , (27) 

where summation over repeated indices is implied. Xij represents the total current-current 
response of graphene to external fields which, decomposed into longitudinal and transverse 
contributions, is given by 

" i=« ~* ^b? t i lJ " • (28) 

Longitudinal and transverse components thus decouple and in the following, we will only 
explicitly label these different channels when it is necessary. The full photon propagator of 
a general layered structure is obtained by superposing all possible scattering paths similar 
to Eq. (27). 

Let us finally recall that the gauge field at any point r due to a current source at r' is 
given through the photon propagator via 

A a (r) = - J d?r'V a P{r,r')f{r') . (29) 

Due to this linear correspondence, the photon propagator can be deduced from the scat- 
tering problem of the gauge field using the expressions of V^f in a homogeneous medium. 

In the following will outline the general scattering problem for one interface, discussing 
both, the longitudinal and transverse channel. We will then give the explicit expressions 
for the retarded photon propagator, dielectric function and graphene loss function for a 
arbitrary double layer structure. 

3.1 Scattering on one interface 

All general properties of the photon propagator of a layered geometry can be deduced 
from the scattering problem of one single interface. We will, therefore, discuss this simple 
problem in some detail. The generalizations are then straightforward. In the following, 
we will set the plane of incidence the xz-plane and the interface at z = 0. 

3.1.1 Longitudinal polarization 

For longitudinal polarization, the general vector field has a component parallel (x) and 
normal (z) to the interface: 

A(r, z) = Y J e iqr (Q, z)e q + A ± (q, z)e z ) . (30) 
q 

With q[ = yjq 2 — (uj/a) 2 and c» the speed of light in the corresponding medium, we make 
the ansatz (j =||, _L) 

*<«•*>-{&■<• :.tl ■ (3i) 
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The components of A 1 - are obtained from the components of A" via the condition for a 
transverse field V • A = 0. This gives the following relations: 

= i^-af , at = -'^4 . at = ^4 (32) 

From Eqs. (24), we see that the parallel component of the vector field is continuous 
at the interface, but the normal component of the displacement field makes a jump if a 
graphene layer is present leading to a charge density p at the interface. This component is 
related to the vector field via the relation D- 1 — esoiojA- 1 . With the continuity equation 
up — q • j = and the linear response j = — XiA q , the set of equations closes. Together 
with Eq. (32), we thus obtain the following two conditions: 



al+al = a\ (33) 
-pXi a t 

The transmission and reflection amplitude for longitudinal polarization then read 



e 2 q[4 - e iq ' 2 (4 ~ 4) = ^kxUl (34) 



T _ 4 _ 2q|ei 



£0^ 



II a' c n't l'll2x"(q,^) 

^=^ = g2£l ' gl£2 + T7^ . (36) 
3.1.2 Transverse polarization 

For transverse polarized light and the plane of incidence again in the icz-plane, only the 
y-componcnt A y of the vector field is non-zero. We thus make the following ansatz: 

From Eqs. (25) we see that the vector potential is continuous at the interface and that 
the first derivative makes a jump due to the current generated by the vector field inside 
the graphene plane. The current is again related to the corresponding transverse current- 
current susceptibility, Xt > y i a linear response. [44, 45] We thus obtain the following two 
conditions: 

a, + a r = a t (38) 
- — at - — (a r - a,i) = poXt a t (39) 

The transmission and reflection amplitude for transverse polarization then read 
rr _ a t _ ^2q[ 



a% M29i +Mi9 2 + MiM2MoX?(<?,^) ' 



(40) 



R= at = M2gj - Mig 2 - MiM2MoX?(<L u) / 41 s 
a, p 2 q[ + Mi<? 2 + MiA i 2A 1 oX?(<?,^) ' 
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Figure 1: (color online): Schematic setup of the double layer graphene structure. The two 
graphene layers, characterized by graphene's current response of the two layers xi and X2 and 
separated by a distance d, are sandwiched by the three dielectric media characterized by e±, 
€2, and £3. The corresponding magnetic permeabilities fi±, fi2, and /J3 are suppressed. 



3.2 Scattering on two interfaces 

We will now discuss in some detail the special case of a double layer graphene structure 
by applying the basic matching conditions discussed in the previous subsection to two 
interfaces. We will give explicit expressions for the retarded photon propagator, the di- 
electric function and the full response of graphene including the definition of the energy 
loss function. 



3.2.1 Photon propagator 

For two and more interfaces it is convenient to write Fourier transformed part of the 
photon propagator as a n x n-matrix, n denoting the number of interfaces. The full 
photon propagator then satisfies the usual Dyson-like equations (one for each polarization 
channel) 

d=(l-d°x°)- 1 d Q , (42) 

where for the special case of two graphene layers the matrix x° = diag(xi)X2) represents 
the bare graphene's response in layer 1 (xf) and layer 2 (x°). dP is the photon propagator 
in the absence of graphene (x° = 0), but with the dielectric geometry of Fig. 1 . 

The entries of the matrix d can be obtained from the standard matching conditions or, 
equivalently, using multiple scattering formalism. In the latter case, they can be written 
as 

(d)n = di (1 + fi, 3 ) (d)ia = d 3 4,1 (43) 

(d)ai = d x t 1)3 (d) 22 = d 3 (I + f 3 ,i), (44) 

where di — 2e .3' oLJ 2 for longitudinal polarization, and di = — ^^P- for the transverse case, 
see Eq. (21). The compound reflection and transmission amplitudes are 



l-rairaae- 2 ^ (4g) 
h 2 t2 3 e- q i d 



113 ~ 1 



S 



with the obvious expression for index exchange, where rij(Uj) are the corresponding coef- 
ficients for a single graphene interface from medium i into medium j, previously obtained. 
Notice that a recursive interpretation of the last formula can be used for calculating the 
reflection and transmission amplitudes for any layered structure with multiple interfaces. 

Solving the matching conditions for the scattering problem with two interfaces but 
without graphene (\ — 0) directly, one can obtain more explicit expressions for the 
photon propagator in the absence of graphene, d°. Using the results of Ref. [32], one 
obtains for longitudinal polarization the following compact result: 

d o = gWk^2 ( cosh(^d) + || sinh(^) 1^ 
1 e co 2 Ni fi y l cosh(^rf) + g sinh(^rf) 

with N L0 = q^W^i + Qi e 3) cosh(q' 2 d) + (g 2 2e i e 3 + qWs^) smh(q 2 d). For the transverse 
part, one obtains the corresponding expression 

d o = _ m^M^A ( cosh ( < l2 d ) + ^ft sinh(g^) 1 
' ' '"" N tfi { 1 " cosh(^d) + ^|sinh(^d) 

with N t ,o = /^(/Wi + cosh(q' 2 d) + {nhWs + MiM3<7 2 2 ) smh(q' 2 d). 

3.2.2 Graphene's response 

Graphene's response obeys similar equations (one for each polarization channel) 

X=(l-X°dT 1 X°, (48) 

which, together with Eq. (42) , provide the complete dynamics of the coupled matter-field 
system. The retarded dielectric function of double layer graphene with nonhomogeneous 
background is then often defined by[31] 

e(q 7 cj)^dct(l-x°d ) . (49) 

Usually — line -1 is used to discuss the plasmonic spectrum, but this function changes sign 
and can thus not be interpreted as (positive definite) spectral density. But instead of the 
determinant, we find it more convenient to discuss the trace of the the full response matrix. 
Graphene's excitations correspond to the imaginary part of the full response, and to reveal 
its presence we will discuss the following generalization of the energy loss function S(q, lj) 
to several layers: 

S(q, u) = -lm X (q, w) = ~lm Tr X (q, w) (50) 

Since S(q, lo) is related to the imaginary part of a causal function, it is strictly positive 
and since it also proportional to the usual definition of the energy loss function for a single 
layer, Eq. (50) can serve as a straightforward generalization of the energy loss function 
for arbitrary layered systems. 

Let us briefly comment on the physical interpretation of the response matrix and the 
related energy loss function. The diagonal entries of the response matrix x are given by 
the response of a particular layer if the gauge field is only applied to just this particular 
layer. Diagonalizing the response matrix \, one can discuss the elementary excitations of 
the full system, separately. This was done in Ref. [35], where for double-layer graphene 
the in-phase and out-of-phase excitations were analyzed. Since the trace of the response 
matrix x is invariant with respect to unitary transformations, it serves as natural choice for 
the definition of the generalized energy loss function and for the discussion of the internal 
excitations of the whole system, i.e., the sum of excitations of all layers. 
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Figure 2: (color online): The energy loss function — Im%(g, a; + iO) in units of ep/h 2 for 
longitudinal polarization for single layer (left) and double layer with kpd = 0.35 (right) at 
a temperature T = Tp/4 with the same dielectric medium for all regions e = 1 (air). Also 
shown the zero-temperature plasmon dispersion (black solid lines). 



As mentioned above, the imaginary part of the general response is related to graphene's 
excitations. In the following section, we will use the developed formalism to discuss the 
spectrum of longitudinal plasmonic excitations for single and double layer structures. For 
this, retardation can formally be neglected. But we stress that the formalism can also 
straightforwardly be used for multiple layer structures as well as to discuss the spectrum 
of transverse plasmonic excitations where retardation effects are crucial. 



4 Plasmons in layered structures at finite temperature 

In this section, we will apply our formalism and discuss the longitudinal response of lay- 
ered structures at finite temperature. To do so, we will use the density-density correlation 
function P°(q,Lu) which is related to the longitudinal component of the current-current 
correlation function X;(<7>k>) via the continuity equation. [45] Within the Dirac-cone ap- 
proximation, this reads 

X$(q,u,) = ^P°(q,u>). (51) 
4.1 Polarizability of graphene 

Up to now, there is no approximation involved except of decomposing the current response 
into a longitudinal and transverse contribution which is well justified for any transitions 
not too close to the van Hove singularity. We will now approximate the full density-density 
correlation function by the non-interacting polarizability[42] 

P o (a lbJ ) _ 9_s9^ f d 2 k fss ' ( , s nF(£ s (fc))-n f (i^'(|k + q|)) 
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Figure 3: (color online): The energy loss function — Im%(g, a; + iO) in units of ep/h 2 for 
longitudinal polarization for single layer (left) and double layer with kpd = 0.35 (right) at a 
temperature T = Tp with the same dielectric medium for all regions e = 1 (air). Also shown 
the zero-temperature plasmon dispersion (black solid lines). 



with E ± (k) = ±hvpk the eigenenergies, tif(E) = (eP^ E ~^> + 1) _1 the Fermi function, 
and g s = g v = 2 the spin and valley degeneracy for graphene. A characteristic difference 
between the polarizability of graphene and that of a two-dimensional electron gas is the 
appearance of the prefactors f ss (k, q) coming from the band-overlap of the wave function 

where ip denotes the angle between k and q. 

At zero temperature, we have [i — tp with ep the Fermi energy and the analytic solution 
of Ref. [42] can be decomposed in several patches corresponding to inter- and intraband 
transitions, respectively. In Figs. 2-6, the three most basic patches are separated by black 
dashed lines. 

At finite temperature the chemical potential is determined with respect to the electronic 
density n by the following relation: 

/oo 
dev{e) [n F (e) - 6(-e)] = n , (54) 
-oo 

with the density-of-states given by v(e) — g s g v \e\/(2irv F ). For an electron-doped system, 
we have e F > fi > 0. At the neutrality point, n = and due to particle-hole symmetry we 
then have fi = 0. For our numerical calculations, we will make use of the semi-analytical 
expression of the polarizability presented in Ref. [43] . 

4.2 Loss function at finite and zero doping 

We will now consider single and double layer structures at finite and zero doping at low 
and high temperatures. To clarify the discussion we will choose a homogeneous medium, 
i.e., ei = e 2 = e 3 = 1. 
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q/k T q/k T 

Figure 4: (color online): The energy loss function — Imx(q,w + iO) in units of ep/fi 2 (er = 
hvpkp) for longitudinal polarization for single layer (left) and double layer with d = 2nm 
(right) at a temperature T = 300K at zero doping n = with the same dielectric medium for 
all regions e = 1 (air). Also shown the zero-temperature plasmon dispersion for finite doping 
with k T = 7^ 2 In 2 (black solid lines). 



We first discuss the plasmon dispersion at finite temperature for suspended single and 
double layer graphene with a layer separation of kpd = 0.35. For an electronic density 
of n = 10 12 cnT 2 , we then have d = 2nm and T F = 1200K, such that T = Tp/4 would 
correspond to approximately room temperature. The energy loss function is shown for 
T = Tp/4 in Fig. 2 and for T = T F for 3. 

We compare the energy loss function with the plasmon dispersion at zero temperature 
by determining dete = 0. In the case of a finite imaginary part of P°(q,uj), i.e., Landau 
damping, we set ImP°=0. This guarantees the convergence of the two plasmonic modes 
for large g.[32] We see that at intermediate temperatures the plasmon dispersion is red- 
shifted compared to the zero temperature solution, whereas for high temperatures it is 
blue-shifted. This follows directly from the behavior of P°(q = 0, u> — > 0) which is related 
to the Drude weight D and defines the plasmon dispersion. [42] At finite temperature, this 
reads 

D = 9b9v fc fJ ( ln(l + e»' k * T ) + ln(l + e~^ kBT )\ (55) 

which is a non-monotonic function of the temperature due to the temperature dependence 
of the chemical potential given in Eq. (54). 

In Fig. 4, we show the plasmon dispersion at the neutrality point n — at T = 300K. 
It was already discussed earlier that also in this limit (damped) plasmon excitations exist 
in single layer graphene due to the weak Landau damping. [46] In fact, there is an analytical 
approximation which uses the zero temperature result of the polarizability[42] with the 
replacement hp — > kx — ^^2 In 2. [47] This agrees with the physical expectation that at 
finite temperature there is a finite electronic density due to thermal broadening of the 
Fermi function. 

The parameter in Fig. 4 corresponds to a thermal wave number kp = 0.061nm -1 and 
electron density n = 1.2 x 10 n cm~ 2 . The energy loss function is given in units of ep/h 2 
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Figure 5: (color online): The energy loss function — Imx(q, uj + iO) in units of ep/h? for longi- 
tudinal polarization of double layer graphene with equal carrier density and layer separation 
k F d = 1.77 at temperature T = 2>/10 (left) and T = T F (right) with e x = 1, e 2 = 6 and 
e3 = 3.8. Also shown the finite-temperature plasmon dispersion obtained by Re(dete) = for 
ImP° = (black lines) and finite ImP° (green lines). 



with €t — hvpkr ~ 36meV. As one can see, the analytical approximation agrees well with 
the maximum of —lmx(q, uj + iO) even in the case of double layer graphene. 

4.3 Loss function for nonhomogeneous dielectric media 

We will now discuss the effect of a nonhomogeneous dielectric media which can lead to 
changes of the plasmon dispersion due to the different photon propagator. In Fig. 5 we 
show the results for a double layer graphene structure for two temperatures T = Tp/10 
(left) and T = Tp (right). We choose SiC>2 with 63 = 3.8 as a substrate and AI2O3 with 
£2 = 6 as a buffer layer. It is further assumed that the top layer is air with e\ = 1. 
The two layers have the same electron density and layer separation kpd = 1.77 which for 
n = 10 12 cm~ 2 corresponds to d = lOnm. 

We use the same parameters as in Ref. [34] and the results of this reference are shown 
as green solid lines, obtained from the real part of dete = 0. Since the authors do not set 
the imaginary part to zero, the in-phase and out-of-phase modes do not merge for q^kp 
in contrary to the physical expectation. The correct way to approximate the plasmon 
dispersion from the zeros of the dielectric function is thus to set the imaginary part of 
the current-current correlation function to zero (black solid lines). Even after including 
temperature broadening and Landau damping, the discussion of Ref. [34] remains incorrect 
since it is based on the response of single layer graphene. 

Apart from the mere (numerical) corrections of the plasmonic dispersion due to a 
nonhomogeneous dielectric background, it is possible to shift spectral weight from the 
out-of-phase to the in-phase mode. One can also change the plasmon dispersion from the 
typical square-root dispersion of charged plasmonic waves to a linear dispersion typical 
for acoustic sound waves. In Figs. 2-4, where the same dielectric constant was chosen 
to be the one of air, the out-of-phase mode is more dominant extending towards larger 
wave numbers. In Fig. 5, we see that the two modes are almost equal in weight merging 
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Figure 6: (color online): The energy loss function — lmx(q,0J + iO) in units of ep/h? and 
temperature T = Tp/4 for longitudinal polarization. Right hand side: Double layer graphene 
with kpd = 0.35 with large dielectric substrate €3 = 300 (ei = 62 = 1). Left hand side: 
topological insulator (g v = g s = 1) with width kpd = 2.37 and €2 = 100 (ei = 1, £3 = 4). Also 
shown the zero-temperature plasmon dispersion (black solid lines). 



in the region of finite Landau damping. If we choose now a substrate with a very large 
dielectric constant, the in-phase mode becomes more dominant whereas the out-of-phase 
mode basically merges with the Dirac cone. Since the latter behavior is also seen in the 
bare graphene response, we can thus state that there is only an in-phase plasmonic mode. 
This is shown on the left hand side of Fig. 6. 

The plasmonic mode is also dominant in the case of a topological insulator where the 
buffer layer is resembled by a strong dielectric medium, i.e., we set e\ = 1, 62 = 100, and 
€3 = 4. This can be seen on the right hand side of Fig. 6 where we use the same parameters 
as in Ref. [33]. They are given by g v = g s = 1 with a sample width of kpd = 2.37 and a 
Fermi velocity of vti = 5 x 10 5 m/s which scales out. 

Let us finally comment on the linear plasmon mode present on the left hand side of 
Fig. 6. In the appendix, we give details of how to derive and approximate the plasmon 
dispersion of the upper layer as u> p — v a q. For this we do not include the full expression 
of the current-current correlation function, but only use the long wavelength limit. This 
yields the simple result for the group velocity of the linear (acoustic) mode 

2a g dkp 

v a = — u - — ^v F , (56) 

£2 

where kp is the Fermi wave number of the upper graphene layer and a g = uc/vf graphene's 
fine-structure constant. This approximation breaks down for small dkp and large £2 
since v a cannot become smaller than the Fermi velocity[32] and a more careful analysis is 
necessary. [33] 

Eq. (56) can nevertheless be used to discuss several aspects. First, v a only depends 
on the buffer substrate; second, only the upper graphene layer enters in the expression 
since the lower one is perfectly screened. This formula can thus not only be applied to 
double layer graphene on substrates with large dielectric constants, but also to single layer 
graphene on top of metals interpreting d as the graphene-metal distance. We can thus 
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calculate the group velocity of the plasmonic mode for a single layer graphene on top 
of Pt(lll) as discussed experimentally in Ref. [38] where the energy loss function was 
measured. With the parameters = 0.3eV, d = 0.33nm and £2 = 1 corresponding to the 
experimental set-up, we have kpd = 0.15 and thus v a = 1.15vp- The screened plasmon 
dispersion will thus lie close to the Dirac "light-cone" in agreement with experiment. 

5 Summary 

We presented results of the energy loss function of layered structures which is the only way 
to unambiguously discuss the plasmon dispersion in the presence of dissipative terms like 
Landau damping or finite temperature. We discussed the effect of temperature and non- 
homogeneous dielectric medium, including large dielectrics which almost perfectly screen 
the graphene layers, but our formalism equally applies to systems with unbalanced elec- 
tronic densities. Our main results are (i) the possibility of shifting relative weight of the 
several plasmonic branches by changing the nonhomogeneous dielectric background and 
(ii) a simple formula for the sound velocity of the (linear) plasmonic mode of the upper 
graphene layer in the case of a substrate with large dielectric constant. These insights 
might be useful towards the engineering of specific plasmon modes for future plasmonic 
circuitries based on graphene. 
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7 Appendix: Large substrate screening and acoustic 
plasmons 

In this appendix, we show that a large (huge) value of e 3 renders medium 3 almost a 
metal, largely screening the Coulomb interaction in the upper graphene layer and turning 
the otherwise square-root plasmon into an acoustic mode over a wide wave number range. 
This is similar to the proposal of Ref. [44] where a perfect metal as substrate is considered. 
The acoustic mode evolves into a regular two-dimensional plasmon only at much reduced 
wave numbers, due to the incomplete screening of medium 3. 

7.1 Linear plasmon mode of the upper graphene layer 

We assume the usual geometry where the three dielectric media separate the upper and 
lower graphene layer, see Fig. 1. We assume a very large value of e 3 s=y 300, much greater 
than the other regular values of ei. 2 , and ignore retardation. Plasmons are solutions of 
(seeEq. (45)) 



1 - r 2 ir 2 3e 



-2qd _ 



= 0. 



(57) 



The reflection amplitude r 2 3 is given by (see Eq. (36) in the limit c — > oo) 



r-23 = 




(58) 



£ 2 + £ 3 - J^J2 
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In the limit £3 » £1,2, ^23 ~ — 1, i-e., there is perfect screening of the lower graphene layer. 
On the other hand r 21 , given by 



r 21 = ^7T, (59) 

£ 2 + £ 1 - 

can be rewritten upon ignoring the q dependence of x? — > TT S ^r LvF ^ F (l° ca l response) as 

e2 — ei , ,2 1 , ,2 

r 2 i = (60) 

- w 2 

where w 2 is the plasmon dispersion for the upper graphene layer between the semi-infinite 
dielectrics £1,2, given by 

^ = 5TTSS* <61) 



Now, the solution of 1 + r 21 e 2qd = in the limit q — > is 

cj 2 -+ £ -^^uj 2 2 q d. (62) 



This is clearly an acoustic mode. The physics is simple: medium 3 acts just as a metal, 
screening the long-range Coulomb interactions in upper graphene layer, turning the original 
square-root plasmon into an acoustic mode. 

7.2 Range of validity: Regular plasmon 

For sufficiently small q, the screening cannot be perfect and the long range nature of 
the interaction should show up anyway. The previous analysis uses the approximation 
r 23 = — 1. One can see that the first correction to lowest order in — is 

r 23 (g = 0,w^0) = -l+ — , (63) 

£3 

and using this result, the range of validity of the acoustic regime can be established as 

- « Qd < 1. (64) 
£3 

For smaller wave numbers qd < j 3 -, one recovers the standard square root behavior due 
to long-range Coulomb interactions, albeit with much reduced frequencies. To the alluded 
order, one easily finds the regular plasmon as 

- 2 «^H for gd<X (65) 
£3 £3 
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